* Revesre beauty simulations with simultaneous noise draws

quietly do "Scholz-Sicinski data prep.do"
drop aa003re-avt06sp sa006re-nb002rer ca003re-dh120 xidof77-xalcm58 ba001re-bvt06sp

foreach var of newlist  pop exp expsq mar tenure union outwis ind ieduc {
qui gen `var'=.
}
foreach var of varlist pop exp expsq mar tenure union outwis ind ieduc {
qui replace `var'=`var'92
}

log using "Scholz-Sicinski reverse simulation corr 025.log", replace

gen simadd=lwage92+0.04*zbeautyf if female==0 & zerowage92==0 & selfemp92==0  //Need to add noise, t-stat too high

quietly xi:regress lwage92 zbeautyf iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0
predict resid92 if e(sample), resid	//obtain residuals
quietly sum resid92,d
scalar lwage92sd=r(sd)			//save the standard deviation

quietly sum lwage92 if female==0 & zerowage92==0 & selfemp92==0,d //obtain the mean of log wages
scalar lwage92mean=r(mean)		//save the mean

*turn below into function
program beausim, rclass


*FIRST STAGE - STRENGTHENING THE RELATIONSHIP BETWEEN EARNINGS AND BEAUTY
gen noise=rnormal(lwage92mean,lwage92sd)
gen simlwage92=simadd+noise if female==0 & zerowage92==0 & selfemp92==0

*SECOND STAGE - ATTENUATION IN BEAUTY
gen b92 = rnormal(0,1)
egen b92n=std(b92)
gen a92=0.25/sqrt(1-0.25^2)
gen b92c=a92*zbeautyf+b92n
egen beau92=std(b92c)

xi:regress simlwage92 beau92 iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0

drop b92 b92n a92 b92c beau92 noise simlwage92

end

simul _b _se ,reps(100000) :beausim

gen t=_b_beau92/_se_beau92

gen onestar=0
replace onestar=1 if t>=1.645 

gen twostar=0
replace twostar=1 if t>=1.96

sum _b_beau92
sum _se_beau92
sum t
tab1 onestar twostar
sum _b_beau92

log c
clear all



quietly do "Scholz-Sicinski data prep.do"
drop aa003re-avt06sp sa006re-nb002rer ca003re-dh120 xidof77-xalcm58 ba001re-bvt06sp

foreach var of newlist  pop exp expsq mar tenure union outwis ind ieduc {
qui gen `var'=.
}
foreach var of varlist pop exp expsq mar tenure union outwis ind ieduc {
qui replace `var'=`var'92
}

log using "Scholz-Sicinski reverse simulation corr 040.log", replace

gen simadd=lwage92+0.04*zbeautyf if female==0 & zerowage92==0 & selfemp92==0  //Need to add noise, t-stat too high

quietly xi:regress lwage92 zbeautyf iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0
predict resid92 if e(sample), resid	//obtain residuals
quietly sum resid92,d
scalar lwage92sd=r(sd)			//save the standard deviation

quietly sum lwage92 if female==0 & zerowage92==0 & selfemp92==0,d //obtain the mean of log wages
scalar lwage92mean=r(mean)		//save the mean

*turn below into function
program beausim, rclass


*FIRST STAGE - STRENGTHENING THE RELATIONSHIP BETWEEN EARNINGS AND BEAUTY
gen noise=rnormal(lwage92mean,lwage92sd)
gen simlwage92=simadd+noise if female==0 & zerowage92==0 & selfemp92==0

*SECOND STAGE - ATTENUATION IN BEAUTY
gen b92 = rnormal(0,1)
egen b92n=std(b92)
gen a92=0.40/sqrt(1-0.40^2)
gen b92c=a92*zbeautyf+b92n
egen beau92=std(b92c)

xi:regress simlwage92 beau92 iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0

drop b92 b92n a92 b92c beau92 noise simlwage92

end

simul _b _se ,reps(100000) :beausim

gen t=_b_beau92/_se_beau92

gen onestar=0
replace onestar=1 if t>=1.645 

gen twostar=0
replace twostar=1 if t>=1.96

sum _b_beau92
sum _se_beau92
sum t
tab1 onestar twostar
sum _b_beau92

log c
clear all


quietly do "Scholz-Sicinski data prep.do"
drop aa003re-avt06sp sa006re-nb002rer ca003re-dh120 xidof77-xalcm58 ba001re-bvt06sp

foreach var of newlist  pop exp expsq mar tenure union outwis ind ieduc {
qui gen `var'=.
}
foreach var of varlist pop exp expsq mar tenure union outwis ind ieduc {
qui replace `var'=`var'92
}

log using "Scholz-Sicinski reverse simulation corr 063.log", replace

gen simadd=lwage92+0.04*zbeautyf if female==0 & zerowage92==0 & selfemp92==0  //Need to add noise, t-stat too high

quietly xi:regress lwage92 zbeautyf iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0
predict resid92 if e(sample), resid	//obtain residuals
quietly sum resid92,d
scalar lwage92sd=r(sd)			//save the standard deviation

quietly sum lwage92 if female==0 & zerowage92==0 & selfemp92==0,d //obtain the mean of log wages
scalar lwage92mean=r(mean)		//save the mean

*turn below into function
program beausim, rclass


*FIRST STAGE - STRENGTHENING THE RELATIONSHIP BETWEEN EARNINGS AND BEAUTY
gen noise=rnormal(lwage92mean,lwage92sd)
gen simlwage92=simadd+noise if female==0 & zerowage92==0 & selfemp92==0

*SECOND STAGE - ATTENUATION IN BEAUTY
gen b92 = rnormal(0,1)
egen b92n=std(b92)
gen a92=0.63/sqrt(1-0.63^2)
gen b92c=a92*zbeautyf+b92n
egen beau92=std(b92c)

xi:regress simlwage92 beau92 iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0

drop b92 b92n a92 b92c beau92 noise simlwage92

end

simul _b _se ,reps(100000) :beausim

gen t=_b_beau92/_se_beau92

gen onestar=0
replace onestar=1 if t>=1.645 

gen twostar=0
replace twostar=1 if t>=1.96

sum _b_beau92
sum _se_beau92
sum t
tab1 onestar twostar
sum _b_beau92

log c
clear all


quietly do "Scholz-Sicinski data prep.do"
drop aa003re-avt06sp sa006re-nb002rer ca003re-dh120 xidof77-xalcm58 ba001re-bvt06sp

foreach var of newlist  pop exp expsq mar tenure union outwis ind ieduc {
qui gen `var'=.
}
foreach var of varlist pop exp expsq mar tenure union outwis ind ieduc {
qui replace `var'=`var'92
}

log using "Scholz-Sicinski reverse simulation corr 087.log", replace

gen simadd=lwage92+0.04*zbeautyf if female==0 & zerowage92==0 & selfemp92==0  //Need to add noise, t-stat too high

quietly xi:regress lwage92 zbeautyf iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0
predict resid92 if e(sample), resid	//obtain residuals
quietly sum resid92,d
scalar lwage92sd=r(sd)			//save the standard deviation

quietly sum lwage92 if female==0 & zerowage92==0 & selfemp92==0,d //obtain the mean of log wages
scalar lwage92mean=r(mean)		//save the mean

*turn below into function
program beausim, rclass


*FIRST STAGE - STRENGTHENING THE RELATIONSHIP BETWEEN EARNINGS AND BEAUTY
gen noise=rnormal(lwage92mean,lwage92sd)
gen simlwage92=simadd+noise if female==0 & zerowage92==0 & selfemp92==0

*SECOND STAGE - ATTENUATION IN BEAUTY
gen b92 = rnormal(0,1)
egen b92n=std(b92)
gen a92=0.87/sqrt(1-0.87^2)
gen b92c=a92*zbeautyf+b92n
egen beau92=std(b92c)

xi:regress simlwage92 beau92 iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0

drop b92 b92n a92 b92c beau92 noise simlwage92

end

simul _b _se ,reps(100000) :beausim

gen t=_b_beau92/_se_beau92

gen onestar=0
replace onestar=1 if t>=1.645 

gen twostar=0
replace twostar=1 if t>=1.96

sum _b_beau92
sum _se_beau92
sum t
tab1 onestar twostar
sum _b_beau92

log c
clear all


* Just first stage 
quietly do "Scholz-Sicinski data prep.do"
drop aa003re-avt06sp sa006re-nb002rer ca003re-dh120 xidof77-xalcm58 ba001re-bvt06sp

foreach var of newlist  pop exp expsq mar tenure union outwis ind ieduc {
qui gen `var'=.
}
foreach var of varlist pop exp expsq mar tenure union outwis ind ieduc {
qui replace `var'=`var'92
}

log using "Scholz-Sicinski reverse simulation corr 100.log", replace

gen simadd=lwage92+0.04*zbeautyf if female==0 & zerowage92==0 & selfemp92==0  //Need to add noise, t-stat too high

quietly xi:regress lwage92 zbeautyf iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0
predict resid92 if e(sample), resid	//obtain residuals
quietly sum resid92,d
scalar lwage92sd=r(sd)			//save the standard deviation

quietly sum lwage92 if female==0 & zerowage92==0 & selfemp92==0,d //obtain the mean of log wages
scalar lwage92mean=r(mean)		//save the mean

*turn below into function
program beausim, rclass


*FIRST STAGE - STRENGTHENING THE RELATIONSHIP BETWEEN EARNINGS AND BEAUTY
gen noise=rnormal(lwage92mean,lwage92sd)
gen simlwage92=simadd+noise if female==0 & zerowage92==0 & selfemp92==0

xi:regress simlwage92 zbeautyf iqscore fathed mothed fathproex pinc57 farmbg i.ipop57 sibsnum  ///
 fullhh momemployed  hssize milty pop exp expsq mar tenure union outwis i.ind i.ieduc ///
 ivocat if female==0 & zerowage92==0 & selfemp92==0

drop noise simlwage92

end

simul _b _se ,reps(100000) :beausim

gen t=_b_zbeautyf/_se_zbeautyf

gen onestar=0
replace onestar=1 if t>=1.645 

gen twostar=0
replace twostar=1 if t>=1.96

sum _b_zbeautyf
sum _se_zbeautyf
sum t
tab1 onestar twostar
sum _b_beau92

log c
clear all
